Background

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com. Sébastien Nusslé and I started this document to collaborate on our data analysis while he was in Switzerland and I was at Berkeley.  

This is a project I have been working on with Kathleen Matthews, Zack Steel, Stephanie Carlson, and Sébastien Nusslé while I was in Stephanie Carlson’s lab https://nature.berkeley.edu/carlsonlab/. Kathleen had gone out every summer for many years to count yellow legged frogs at Dusy Basin in Kings Canyon National Park (Rana sierra, the Sierra Nevada yellow-legged frog). We compiled counts from 1997-2012. Volunteers went out at least three times a year (July, August, September) and counted the frogs at their various life-stages. We analyzed this dataset as a collaborative effort.

Acknowledgments

First, we would like to thank the summer field assistants (Adam Calo, Alex Figueroa, Ambrose Tuscano, Annie Shattuck, Carrie Sendak, Kyung Koh, Blake Meneken, Christopher Loomis, David Court, Ed Stauber, Elizabeth Gruenstein, Eva Patton, Igor Lacan, Jacob Martin, Jacquelyn Kracke, Jeff Schlueter, Jill Andersen, Joe Fontaine, Julia Anderson, Karen Pope, Krishna Feldman, Leah Culp, Linda Vance, Matthew Young, Michael Davis, Nick Earley, Rebecca Doubledee, Rick Ziegler, Robert Grasso, Steve Hasslinger, and Vincent King who did the surveys from 1997-2012. We would also like to thank Prof. Steve Beissinger (ORCID: 0000-0003-1323-2727) for his advice on the abundance analyses. Timothy Szeliga from NOAA’s Geo-Intelligence Division, William Brown from NOAA’s National Center for Environmental Information, and David Rizzardo from the California Department of Water Resources were especially helpful in finding publicly available precipitation data online. We are also thankful to Cassandra L. Ettinger (ORCID: 0000-0001-7334-403X), Sonia L. Ghose (ORCID: 0000-0001-5667-6876), and Marina E. De León (ORCID: 0000-0001-9973-1951) for their comments on earlier versions of this manuscript and the photograph of a Sierra Nevada yellow-legged frog in Fig. 1.

 

This script assumes you have a working directory in your home directory. Here, this directory is called “/DusyBasin”

All the input files are in “~/DusyBasin/Input_summaryStats/” ADD URL GITHUB.

 

Set-up R and all the necessary libraries:

Take hash away if you still need to install these libraries!

# Take hash away if you still need to install these libraries!

#install.packages(corrplot)
#install.packages(scales)
#install.packages(lme4)
#install.packages(lmerTest)
#install.packages(unmarked)

library(corrplot)
library(scales)
library(lme4)
library(lmerTest)
library(unmarked)

p.trans <- function(p) {
    res <- as.numeric(NA, length(p))
    for (i in 1:length(p)) {
        if (p[i] > 0.1) res[i] <- "p > 0.05 [NS]"
        if (p[i] > 0.05 & p[i] <= 0.1) res[i] <- "p < 0.1 [.]"
        if (p[i] > 0.01 & p[i] <= 0.05) res[i] <- "p < 0.05 [*]"
        if (p[i] > 0.001 & p[i] <= 0.01) res[i] <- "p < 0.01 [**]"
        if (p[i] <= 0.001) res[i] <- "p < 0.001 [***]"
    }
    return(res)}

# Data
load(file = "~/DusyBasin/Input_summaryStats/data.frogs.Rdata")

First, we need a nice map for Figure 1!

Code for Figure 1:

Some more code to play with maps. It is not executed here. Change eval=TRUE, echo=TRUE and it will run.

Now let’s dive into the analysis!

Data exploration…

# Selection of meaningful lakes only (30 and above are rivers)
data <- subset(data, site %in% c(1:11, 13:29))

# Defining new categories
data$frogs <- data$subadults + data$adults
data$juveniles <- data$metamorphs + data$premetamorphs
data$tadpoles <- data$larvae + data$juveniles

plot(log(data$eggs+1) ~ log(data$larvae+1))
pred <- loess(log(data$eggs+1) ~ log(data$larvae+1))
points(pred$fitted[order(data$larvae)] ~ log(data$larvae+1)[order(data$larvae)], type = "l", lwd = 2, col = "red")

# New dataset: pooling seasons per year and sites as mean numbers

eg_mean_y <- tapply(data$eggs, list(data$year, data$site), mean, na.rm = TRUE)
la_mean_y <- tapply(data$larvae, list(data$year, data$site), mean, na.rm = TRUE)
ju_mean_y <- tapply(data$juveniles, list(data$year, data$site), mean, na.rm = TRUE)
su_mean_y <- tapply(data$subadults, list(data$year, data$site), mean, na.rm = TRUE)
ad_mean_y <- tapply(data$adults, list(data$year, data$site), mean, na.rm = TRUE)
fr_mean_y <- tapply(data$frogs, list(data$year, data$site), mean, na.rm = TRUE)
ta_mean_y <- tapply(data$tadpoles, list(data$year, data$site), mean, na.rm = TRUE)

ndat <- data.frame(year = as.numeric(rep(rownames(eg_mean_y), times = dim(eg_mean_y)[2])), site = as.numeric(rep(colnames(eg_mean_y), each = dim(eg_mean_y)[1])), eggs = as.numeric(eg_mean_y), larvae = as.numeric(la_mean_y), tadpoles = as.numeric(ta_mean_y), juveniles = as.numeric(ju_mean_y), subadults = as.numeric(su_mean_y), adults = as.numeric(ad_mean_y), frogs = as.numeric(fr_mean_y))
summary(ndat)
##       year           site             eggs             larvae       
##  Min.   :1997   Min.   : 1.000   Min.   : 0.0000   Min.   :   0.00  
##  1st Qu.:2000   1st Qu.: 4.000   1st Qu.: 0.0000   1st Qu.:   0.00  
##  Median :2005   Median : 7.500   Median : 0.0000   Median :  11.71  
##  Mean   :2005   Mean   : 9.214   Mean   : 0.8973   Mean   : 114.82  
##  3rd Qu.:2009   3rd Qu.:11.000   3rd Qu.: 0.1250   3rd Qu.: 100.20  
##  Max.   :2012   Max.   :22.000   Max.   :25.0000   Max.   :2045.80  
##                                  NA's   :19        NA's   :19       
##     tadpoles         juveniles        subadults           adults       
##  Min.   :   0.00   Min.   :  0.00   Min.   :  0.000   Min.   :  0.000  
##  1st Qu.:   0.00   1st Qu.:  0.00   1st Qu.:  0.000   1st Qu.:  0.000  
##  Median :  11.71   Median :  0.00   Median :  0.800   Median :  2.250  
##  Mean   : 137.95   Mean   : 22.96   Mean   :  6.748   Mean   :  8.656  
##  3rd Qu.: 131.00   3rd Qu.:  4.75   3rd Qu.:  4.589   3rd Qu.:  8.589  
##  Max.   :2046.20   Max.   :427.89   Max.   :115.333   Max.   :101.000  
##  NA's   :19        NA's   :19       NA's   :19        NA's   :19       
##      frogs        
##  Min.   :  0.000  
##  1st Qu.:  0.450  
##  Median :  3.333  
##  Mean   : 15.446  
##  3rd Qu.: 17.278  
##  Max.   :170.000  
##  NA's   :19

Code for Figure 2

Let’s move on to proportion of occupied lakes:

# Year goodness of fit: proportion of lakes with eggs, larvae, juveniles, or adults

YGFeggs <- eg_mean_y
    YGFeggs[eg_mean_y > 0] <- 1
    YGFeggs <- apply(YGFeggs, 1, mean, na.rm = TRUE)
LGFeggs <- eg_mean_y
    LGFeggs[eg_mean_y > 0] <- 1
    LGFeggs <- apply(LGFeggs, 2, mean, na.rm = TRUE)

YGFlarvae <- la_mean_y
    YGFlarvae[la_mean_y > 0] <- 1
    YGFlarvae <- apply(YGFlarvae, 1, mean, na.rm = TRUE)

YGFjuveniles <- la_mean_y
    YGFjuveniles[la_mean_y > 0] <- 1
    YGFjuveniles <- apply(YGFjuveniles, 1, mean, na.rm = TRUE)

YGFadults <- ad_mean_y
    YGFadults[ad_mean_y > 0] <- 1
    YGFadults <- apply(YGFadults, 1, mean, na.rm = TRUE)

    
YGF <- rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults)
rownames(YGF) <- c("Eggs", "Tadpoles", "Subadults", "Adults")


# First I have to add three more years with no counts to the plotting dataset:

more.years <- matrix(c(0,0,0,0,0,0,0,0,0,0,0,0), nrow=4, ncol=3, byrow=TRUE)
rownames(more.years) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
colnames(more.years) <- c("2013", "2014", "2015")

YGF <- as.matrix(YGF)
test <- cbind(YGF, more.years)

YGF <- test

par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF, beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = terrain.colors(4), ylab = "Proportion of occupied lakes")

I also did an additional plot where I only included egg masses that were counted during the breeding season!

data.breeding <- subset(data, data$season=="b")
eg_mean_y_b <- tapply(data.breeding$eggs, list(data.breeding$year, data.breeding$site), mean, na.rm = TRUE)

YGFeggs <- eg_mean_y_b
    YGFeggs[eg_mean_y_b > 0] <- 1
    YGFeggs <- apply(YGFeggs, 1, mean, na.rm = TRUE)

YGFlarvae <- la_mean_y
    YGFlarvae[la_mean_y > 0] <- 1
    YGFlarvae <- apply(YGFlarvae, 1, mean, na.rm = TRUE)

YGFjuveniles <- la_mean_y
    YGFjuveniles[la_mean_y > 0] <- 1
    YGFjuveniles <- apply(YGFjuveniles, 1, mean, na.rm = TRUE)

YGFadults <- ad_mean_y
    YGFadults[ad_mean_y > 0] <- 1
    YGFadults <- apply(YGFadults, 1, mean, na.rm = TRUE)

    
YGF <- rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults)
## Warning in rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults): number of
## columns of result is not a multiple of vector length (arg 1)
rownames(YGF) <- c("Eggs", "Tadpoles", "Subadults", "Adults")


# First I have to add three more years with no counts to the plotting dataset:

more.years <- matrix(c(0,0,0,0,0,0,0,0,0,0,0,0), nrow=4, ncol=3, byrow=TRUE)
rownames(more.years) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
colnames(more.years) <- c("2013", "2014", "2015")

YGF <- as.matrix(YGF)
test <- cbind(YGF, more.years)

YGF <- test

par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF, beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = terrain.colors(4), ylab = "Proportion of occupied lakes")

And one with different colors:

YGFeggs <- eg_mean_y
    YGFeggs[eg_mean_y > 0] <- 1
    YGFeggs <- apply(YGFeggs, 1, mean, na.rm = TRUE)
LGFeggs <- eg_mean_y
    LGFeggs[eg_mean_y > 0] <- 1
    LGFeggs <- apply(LGFeggs, 2, mean, na.rm = TRUE)

YGFlarvae <- la_mean_y
    YGFlarvae[la_mean_y > 0] <- 1
    YGFlarvae <- apply(YGFlarvae, 1, mean, na.rm = TRUE)

YGFjuveniles <- la_mean_y
    YGFjuveniles[la_mean_y > 0] <- 1
    YGFjuveniles <- apply(YGFjuveniles, 1, mean, na.rm = TRUE)

YGFadults <- ad_mean_y
    YGFadults[ad_mean_y > 0] <- 1
    YGFadults <- apply(YGFadults, 1, mean, na.rm = TRUE)

    
YGF <- rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults)
rownames(YGF) <- c("Eggs", "Tadpoles", "Subadults", "Adults")


# First I have to add three more years with no counts to the plotting dataset:

more.years <- matrix(c(0,0,0,0,0,0,0,0,0,0,0,0), nrow=4, ncol=3, byrow=TRUE)
rownames(more.years) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
colnames(more.years) <- c("2013", "2014", "2015")

YGF <- as.matrix(YGF)
test <- cbind(YGF, more.years)

YGF <- test

par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF, beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = c("white", "green", "cornflowerblue", "coral1"), ylab = "Proportion of occupied lakes")

Now let’s do a combination of the two plots where I take the sums for egg mass but otherwise the data from the first figure with averages over the three sampling efforts each summer:

# for eggs take only breeding season:
data.breeding <- subset(data, data$season=="b")
eg_mean_y_b <- tapply(data.breeding$eggs, list(data.breeding$year, data.breeding$site), mean, na.rm = TRUE)

YGFeggs <- eg_mean_y_b
    YGFeggs[eg_mean_y_b > 0] <- 1
    YGFeggs <- apply(YGFeggs, 1, mean, na.rm = TRUE)

YGFlarvae <- la_mean_y
    YGFlarvae[la_mean_y > 0] <- 1
    YGFlarvae <- apply(YGFlarvae, 1, mean, na.rm = TRUE)

YGFjuveniles <- la_mean_y
    YGFjuveniles[la_mean_y > 0] <- 1
    YGFjuveniles <- apply(YGFjuveniles, 1, mean, na.rm = TRUE)

YGFadults <- ad_mean_y
    YGFadults[ad_mean_y > 0] <- 1
    YGFadults <- apply(YGFadults, 1, mean, na.rm = TRUE)

    
YGF <- rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults)
## Warning in rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults): number of
## columns of result is not a multiple of vector length (arg 1)
rownames(YGF) <- c("Eggs", "Tadpoles", "Subadults", "Adults")


# First I have to add three more years with no counts to the plotting dataset:

more.years <- matrix(c(0,0,0,0,0,0,0,0,0,0,0,0), nrow=4, ncol=3, byrow=TRUE)
rownames(more.years) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
colnames(more.years) <- c("2013", "2014", "2015")

YGF <- as.matrix(YGF)
test <- cbind(YGF, more.years)

YGF <- test

YGF[1,2]<- 0.7
YGF[1,13]<- 0
YGF[1,14]<- 0
YGF[1,15]<- 0

par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF, beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = terrain.colors(4), ylab = "Proportion of occupied lakes")

Figures 3 and 4 were done using code in another script (directory Scripts_unmarked on GitHub repository https://github.com/megaptera-helvetiae/DusyBasin) by Zack Steel!!!

ADD URL GITHUB.

Snow Data

snow2 <- subset(season_length, season == "winter")[,-4]
    names(snow2) <- c("w.start", "w.end", "year", "w.dur")
    snow2$year <- 1989:2016
snow3 <- subset(season_length, season == "summer")[,-4]; names(snow3) <- c("s.start", "s.end", "year", "s.dur")
    names(snow3) <- c("s.start", "s.end", "year", "s.dur")
    snow3$year <- 1989:2017
snow2 <- merge(snow2, snow3, all = TRUE)

snow2$snow.max <- tapply(dsnow$watcont, dsnow$year, max, na.rm = TRUE)[-1]

par(mfrow = c(2,2), las = 1, mar = c(3,4,1,1))
plot(watcont ~ date, dsnow, type = "l", ylab = "Water Content Equivalent")
plot(tapply(dsnow$watcont, dsnow$year, max, na.rm = TRUE) ~ c(1988:2017), type = "b", pch = 17, ylab = "Maximal snow")
plot(duration ~ year, subset(season_length, season == "winter"), type = "b", pch = 16, ylab = "Winter duration")
    abline(lm(duration ~ year, subset(season_length, season == "winter")))
plot(duration ~ year, subset(season_length, season == "summer"), type = "b", pch = 17, ylab = "Summer duration")
    abline(lm(duration ~ year, subset(season_length, season == "summer")))

Code for Figure 5

Snow plots

par(las = 1, mar = c(5,4,1,1), mfrow = c(2,1))
plot(SnHe_cm ~ Year, snow, type = "b", pch = 16)
boxplot(SnHe_cm ~ Year, snow, type = "b", pch = 16)

# Snow resampling (1 value per year)

boot = 100
years <- unique(snow$Year)
bootsnow <- matrix(NA, boot, length(years))
delta_snow <- NA

t0 <- Sys.time()
for (i in 1:boot) {
    for (j in 1:length(years)) {
        interm <- subset(snow, Year == years[j])$SnHe_cm
        if (length(interm) > 1) bootsnow[i,j] <- sample(x = as.numeric(subset(snow, Year == years[j])$SnHe_cm), size = 1)
        else bootsnow[i,j] <- interm
    }
    res <- lm(bootsnow[i,] ~ years)
    delta_snow[i] <- coef(res)[2]
    }
Sys.time()-t0
## Time difference of 1.214326 secs
hist(delta_snow)
mean(delta_snow); sd(delta_snow)
## [1] -0.5816061
## [1] 0.1461931
t.test(delta_snow)
## 
##  One Sample t-test
## 
## data:  delta_snow
## t = -39.783, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6106139 -0.5525982
## sample estimates:
##  mean of x 
## -0.5816061
# Snow averaging (1 value per year)

snow.year <- tapply(snow$SnHe_cm, snow$Year, mean)
    snow.year2 <- tapply(snow$SnHe_cm, snow$year2, mean)
snow.min <- tapply(snow$SnHe_cm, snow$Year, min)
    snow.min2 <- tapply(snow$SnHe_cm, snow$year2, min)
snow.max <- tapply(snow$SnHe_cm, snow$Year, max)
    snow.max2 <- tapply(snow$SnHe_cm, snow$year2, max)
    
par(las = 1, mar = c(4,4,1,1))
plot(snow.year ~ as.numeric(names(snow.year)), pch = 16, type = "n", ylim = c(0, 500), ylab = "Snow height [cm]", xlab = "Time [year]")
abline(v = c(194:202)*10, lty = 2)
polygon(c(years, rev(years)), c(snow.min, rev(snow.max)), col = "grey")
points(snow.year ~ as.numeric(names(snow.year)), pch = 16, type = "b")
abline(h = 100)
arrows(1961.5, 80, 1974.5, 80, lwd = 2, code = 3)
arrows(1977.5, 80, 1986.5, 80, lwd = 2, code = 3)
arrows(1992.5, 80, 2000.5, 80, lwd = 2, code = 3)
arrows(2001.5, 80, 2006.5, 80, lwd = 2, code = 3)
arrows(2007.5, 80, 2011.5, 80, lwd = 2, code = 3)
text(c(1967.5, 1981.5, 1996, 2004, 2010), y = 70, c(13,9,8,5, 4))

snowheight <- data.frame(years, snow.year, snow.min, snow.max)
res <- lm(snow.year ~ years, subset(snowheight, years %in% 1958:2015));summary(res)
## 
## Call:
## lm(formula = snow.year ~ years, data = subset(snowheight, years %in% 
##     1958:2015))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -139.26  -62.46  -15.95   50.85  198.11 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1655.6963  1276.5364   1.297    0.200
## years         -0.7394     0.6426  -1.151    0.255
## 
## Residual standard error: 81.92 on 56 degrees of freedom
## Multiple R-squared:  0.0231, Adjusted R-squared:  0.005654 
## F-statistic: 1.324 on 1 and 56 DF,  p-value: 0.2547
res <- lm(snow.min ~ years, subset(snowheight, years %in% 1958:2015));summary(res)
## 
## Call:
## lm(formula = snow.min ~ years, data = subset(snowheight, years %in% 
##     1958:2015))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -114.814  -56.987   -4.765   48.538  203.331 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1153.7636  1100.5153   1.048    0.299
## years         -0.5074     0.5540  -0.916    0.364
## 
## Residual standard error: 70.63 on 56 degrees of freedom
## Multiple R-squared:  0.01476,    Adjusted R-squared:  -0.002834 
## F-statistic: 0.8389 on 1 and 56 DF,  p-value: 0.3636
res <- lm(snow.max ~ years, subset(snowheight, years %in% 1958:2015));summary(res)
## 
## Call:
## lm(formula = snow.max ~ years, data = subset(snowheight, years %in% 
##     1958:2015))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -157.59  -75.55  -21.54   81.39  239.97 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2002.1621  1506.2239   1.329    0.189
## years         -0.8956     0.7582  -1.181    0.242
## 
## Residual standard error: 96.67 on 56 degrees of freedom
## Multiple R-squared:  0.02431,    Adjusted R-squared:  0.006889 
## F-statistic: 1.395 on 1 and 56 DF,  p-value: 0.2425
res <- lm(snow.year2 ~ as.numeric(names(snow.year2)));summary(res)
## 
## Call:
## lm(formula = snow.year2 ~ as.numeric(names(snow.year2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -73.070 -22.493   0.649  22.771  55.006 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                   1644.0277   757.5098   2.170   0.0464 *
## as.numeric(names(snow.year2))   -0.7301     0.3835  -1.904   0.0763 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.73 on 15 degrees of freedom
## Multiple R-squared:  0.1946, Adjusted R-squared:  0.1409 
## F-statistic: 3.624 on 1 and 15 DF,  p-value: 0.07632
res <- lm(snow.min2 ~ as.numeric(names(snow.year2)));summary(res)
## 
## Call:
## lm(formula = snow.min2 ~ as.numeric(names(snow.year2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -68.524 -15.743   2.045  19.927 106.525 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   3829.7157   803.7712   4.765 0.000251 ***
## as.numeric(names(snow.year2))   -1.8863     0.4069  -4.635 0.000324 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 41.1 on 15 degrees of freedom
## Multiple R-squared:  0.5889, Adjusted R-squared:  0.5615 
## F-statistic: 21.49 on 1 and 15 DF,  p-value: 0.0003236
res <- lm(snow.max2 ~ as.numeric(names(snow.year2)));summary(res)
## 
## Call:
## lm(formula = snow.max2 ~ as.numeric(names(snow.year2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -157.22  -32.08  -10.32   39.00  147.94 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)                   -1585.8358  1344.8631  -1.179    0.257
## as.numeric(names(snow.year2))     0.9728     0.6809   1.429    0.174
## 
## Residual standard error: 68.77 on 15 degrees of freedom
## Multiple R-squared:  0.1198, Adjusted R-squared:  0.06111 
## F-statistic: 2.041 on 1 and 15 DF,  p-value: 0.1736
plot(snow.year2 ~ as.numeric(names(snow.year2)), pch = 16, type = "b", ylim = c(0, 500))
points(snow.min2 ~ as.numeric(names(snow.year2)), type = "l")
points(snow.max2 ~ as.numeric(names(snow.year2)), type = "l")

#-----------------------------------
# Graphs
#-----------------------------------

boxplot(SnHe_cm ~ Year, subset(snow, Year > 1900))
abline(h = 100, col = "red")

#quartz("", 10, 7)
par(mfrow = c(1,1), las = 1, mar = c(2,4,0,1), oma = c(1,1,1,1))
plot(SnHe_cm ~ Time, snow, type = "n", ylab = "Snow cover [cm]", pch = 16)
abline(h = c(1:10)*50, lty = 2, col = "grey")
abline(v = c(194:202)*10, lty = 2, col = "grey")

points(SnHe_cm ~ Time, snow, type = "p", ylab = "Snow cover [cm]", pch = 16)
points(lowess(snow$SnHe_cm)$y ~ snow$Time, col = "red", lwd = 2, type = "l")
res <- lm(SnHe_cm ~ Time, snow); summary(res)
## 
## Call:
## lm(formula = SnHe_cm ~ Time, data = snow)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -164.12  -65.24  -17.59   60.73  273.98 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 1330.3285   637.8121   2.086   0.0382 *
## Time          -0.5717     0.3216  -1.778   0.0769 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 92.19 on 209 degrees of freedom
## Multiple R-squared:  0.0149, Adjusted R-squared:  0.01018 
## F-statistic: 3.161 on 1 and 209 DF,  p-value: 0.07689
abline(res)
legend("topleft", p.trans(coef(summary(res))[2,4]), bty = "n")

par(mfrow = c(3,1), las = 1, mar = c(2,4,0,1), oma = c(1,1,1,1))
plot(snow.year ~ as.numeric(names(snow.year)), type = "b", ylab = "Mean snow cover")
points(lowess(snow.year)$y ~ as.numeric(names(snow.year)), col = "red", lwd = 2, type = "l")
res <- lm(snow.year ~ as.numeric(names(snow.year))); summary(res)
## 
## Call:
## lm(formula = snow.year ~ as.numeric(names(snow.year)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -140.93  -58.35  -11.23   43.64  197.29 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                  1359.0107   722.4623   1.881   0.0636 .
## as.numeric(names(snow.year))   -0.5885     0.3654  -1.611   0.1113  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76.89 on 79 degrees of freedom
## Multiple R-squared:  0.03179,    Adjusted R-squared:  0.01953 
## F-statistic: 2.594 on 1 and 79 DF,  p-value: 0.1113
abline(res)
legend("topleft", p.trans(coef(summary(res))[2,4]), bty = "n")

plot(snow.min ~ as.numeric(names(snow.year)), type = "b", ylab = "Min snow cover")
points(lowess(snow.min)$y ~ as.numeric(names(snow.year)), col = "red", lwd = 2, type = "l")
res <- lm(snow.min ~ as.numeric(names(snow.year))); summary(res)
## 
## Call:
## lm(formula = snow.min ~ as.numeric(names(snow.year)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -129.419  -51.061   -3.118   44.789  192.067 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  2269.2859   656.1875   3.458 0.000879 ***
## as.numeric(names(snow.year))   -1.0643     0.3319  -3.207 0.001939 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 69.84 on 79 degrees of freedom
## Multiple R-squared:  0.1152, Adjusted R-squared:  0.104 
## F-statistic: 10.28 on 1 and 79 DF,  p-value: 0.001939
abline(res)
legend("topleft", p.trans(coef(summary(res))[2,4]), bty = "n")

plot(snow.max ~ as.numeric(names(snow.year)), type = "b", ylab = "Max snow cover")
points(lowess(snow.max)$y ~ as.numeric(names(snow.year)), col = "red", lwd = 2, type = "l")
res <- lm(snow.max ~ as.numeric(names(snow.year))); summary(res)
## 
## Call:
## lm(formula = snow.max ~ as.numeric(names(snow.year)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -146.41  -63.85  -22.69   58.52  255.04 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                  482.8551   843.0128   0.573    0.568
## as.numeric(names(snow.year))  -0.1317     0.4264  -0.309    0.758
## 
## Residual standard error: 89.72 on 79 degrees of freedom
## Multiple R-squared:  0.001206,   Adjusted R-squared:  -0.01144 
## F-statistic: 0.09539 on 1 and 79 DF,  p-value: 0.7582
abline(res)
legend("topleft", p.trans(coef(summary(res))[2,4]), bty = "n")

Now my favorite plot; Fig. S1 in the manuscript

The summary plot using bubbles to show abundance of different life stages each each. Sites are placed according to their geographic location in the plot.

d <- as.dist(dist)
mds.coor <- cmdscale(d)
mds.coor <- mds.coor[c(1,2,4,5,6,7,9,13),]

par(mar = c(0,0,0,0), mfrow = c(3,5))
for (i in 1:dim(ad_mean_y)[1]) {
    plot(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = la_mean_y[i,c(1,2,4,5:7,9,12)]/50, bg = alpha("green", 0.5), xlim = c(-350, 400), ylim = c(-350, 250), xaxt = "n", yaxt = "n")
    points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ju_mean_y[i,c(1,2,4,5:7,9,12)]/25, bg = alpha("blue", 0.5))
    points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ad_mean_y[i,c(1,2,4,5:7,9,12)]/5, bg = alpha("red", 0.5))
    text(mds.coor[,1], -mds.coor[,2], rownames(mds.coor), font = 2)
    abline(h=0,v=0, col = "grey", lty = 2)
    legend("topright", rownames(la_mean_y)[i], bty = "n", text.font = 4)
}

    #legend("bottomright", legend = c("Larvae (x10)", "Juveniles (x5)", "Adults"), fill = alpha(c("green", "blue", "red"), 0.5), bty = "n")


# And with a legend:

d <- as.dist(dist)
mds.coor <- cmdscale(d)
mds.coor <- mds.coor[c(1,2,4,5,6,7,9,13),]

#quartz("", 11, 7)
par(mar = c(0,0,0,0), mfrow = c(3,5))
for (i in 1:dim(ad_mean_y)[1]) {
  plot(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = la_mean_y[i,c(1,2,4,5:7,9,12)]/50, bg = alpha("green", 0.5), xlim = c(-350, 400), ylim = c(-350, 210), xaxt = "n", yaxt = "n")
  points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ju_mean_y[i,c(1,2,4,5:7,9,12)]/21, bg = alpha("blue", 0.5))
  points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ad_mean_y[i,c(1,2,4,5:7,9,12)]/5, bg = alpha("red", 0.5))
  a <- mds.coor[,1]+25
  b <- -mds.coor[,2]+25
  text(a,b, rownames(mds.coor), font = 2, cex=1.2)
  abline(h=0,v=0, col = "grey", lty = 2)
  legend("topright", rownames(la_mean_y)[i], bty = "n", text.font = 4)
  legend("bottomright", legend = c("Tadpoles (x10)", "Subadults (x5)", "Adults", "Tadpoles (x10) & Adults", "Subadults (x5) & Adults", "Everything"), fill = alpha(c("green", "blue", "red", "saddlebrown", "mediumslateblue", "maroon4"), 0.5), bty = "n")
}

# Can I place the legend somewhere else?

d <- as.dist(dist)
mds.coor <- cmdscale(d)
mds.coor <- mds.coor[c(1,2,4,5,6,7,9,13),]

par(mar = c(0,0,0,0), mfrow = c(3,5))


for (i in 1:dim(ad_mean_y)[1]) {
  plot(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = la_mean_y[i,c(1,2,4,5:7,9,12)]/50, bg = alpha("green", 0.5), xlim = c(-350, 400), ylim = c(-350, 210), xaxt = "n", yaxt = "n")
  points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ju_mean_y[i,c(1,2,4,5:7,9,12)]/21, bg = alpha("blue", 0.5))
  points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ad_mean_y[i,c(1,2,4,5:7,9,12)]/5, bg = alpha("red", 0.5))
  a <- mds.coor[,1]+25
  b <- -mds.coor[,2]+25
  text(a,b, rownames(mds.coor), font = 2, cex=1.2)
  abline(h=0,v=0, col = "grey", lty = 2)
  legend("topright", inset=c(0,0), rownames(la_mean_y)[i], bty = "n", text.font = 4)
}
legend("bottomright", inset=c(0,0), legend = c("Tadpoles (x10)", "Subadults (x5)", "Adults", "T & A", "S & A", "All"), fill = alpha(c("green", "blue", "red", "saddlebrown", "mediumslateblue", "maroon4"), 0.5), bty = "n")

# A little bit more polishing:

d <- as.dist(dist)
mds.coor <- cmdscale(d)
mds.coor <- mds.coor[c(1,2,4,5,6,7,9,13),]

par(mar = c(0,0,0,0), mfrow = c(5,3))


for (i in 1:dim(ad_mean_y)[1]) {
  plot(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = la_mean_y[i,c(1,2,4,5:7,9,12)]/50, bg = alpha("green", 0.5), xlim = c(-350, 450), ylim = c(-400, 270), xaxt = "n", yaxt = "n")
  points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ju_mean_y[i,c(1,2,4,5:7,9,12)]/21, bg = alpha("blue", 0.5))
  points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ad_mean_y[i,c(1,2,4,5:7,9,12)]/5, bg = alpha("red", 0.5))
  a <- mds.coor[,1]+25
  b <- -mds.coor[,2]+25
  text(a,b, rownames(mds.coor), font = 2, cex=1.2)
  abline(h=0,v=0, col = "grey", lty = 2)
  legend("topright", inset=c(0,0), rownames(la_mean_y)[i], bty = "n", text.font = 4, cex = 1.5)
}
legend("bottomright", inset=c(0,0), legend = c("Tadpoles (x10)", "Subadults (x5)", "Adults", "T & A", "S & A", "All"), fill = alpha(c("green", "blue", "red", "saddlebrown", "mediumslateblue", "maroon4"), 0.5), bty = "n")

If I plot this last plot in R-Studio it looks quite OK. Legend only in the last plot and years clearly visible. In R Markdown it looks not so nice.  

Now let’s make some useful summary plots. We decided to concentrate on the following life-stages: number of egg masses, tadpoles, subadults and adults.

 

I will plot tadpoles, subadults and adults in one plot and then egg masses in a separate plot. For tadpoles, subadults and adults I also show the log-transformed counts. For egg masses I don’t because there are so many zeros and the log of zero is an undefined number! –> Also: for eggs I only take counts during the breeding season!

Code for Supplementary Figure S4

# plotting with log

par(mfrow = c(1,3), oma=c(0,0,2,0))
boxplot(larvae ~ year, subset(ndat, larvae > 0), log = "y", ylab="", main = "Tadpoles (when present)", las = 2, cex.axis=0.75)
mtext("LOG(Mean counts per year)", side=2, line=3, cex=0.8)
boxplot(subadults ~ year, subset(ndat, subadults > 0), log = "y", main = "Subadults (when present)", las = 2)
boxplot(adults ~ year, subset(ndat, adults > 0), log = "y", main = "Adults (when present)", las = 2)

# trying to plot all years, without LOG, real counts, until 2015:

write.table(ndat, file="ndat_raw.txt")
# --> I manually added three years with zero counts in Excel!
ndat.v2 <- read.delim("~/DusyBasin/Input_summaryStats/ndat_2015.txt")


par(mfrow = c(1,3), oma=c(0,0,2,0))
boxplot(larvae ~ year, data=ndat.v2, las = 2, ylab="", main="Tadpoles", cex.axis=0.75)
mtext("Mean counts per year", side=2, line=3, cex=0.8)
boxplot(subadults ~ year, data=ndat.v2, main = "Subadults", las = 2)
boxplot(adults ~ year, data=ndat.v2, main = "Adults", las = 2)

Code for Supplementary Figures S5 and S7

I am plotting the eggs separately.

# for eggs take only breeding season:
data.breeding <- subset(ndat.v2, data$season=="b")
eg_mean_y_b <- tapply(data.breeding$eggs, list(data.breeding$year, data.breeding$site), mean, na.rm = TRUE)

ndat.egg <- data.frame(year = as.numeric(rep(rownames(eg_mean_y_b), times = dim(eg_mean_y_b)[2])), site = as.numeric(rep(colnames(eg_mean_y_b), each = dim(eg_mean_y_b)[1])), eggs = as.numeric(eg_mean_y_b))
summary(ndat.egg)
##       year           site             eggs        
##  Min.   :2002   Min.   : 1.000   Min.   : 0.0000  
##  1st Qu.:2003   1st Qu.: 4.000   1st Qu.: 0.0000  
##  Median :2006   Median : 7.500   Median : 0.0000  
##  Mean   :2007   Mean   : 9.214   Mean   : 0.9900  
##  3rd Qu.:2013   3rd Qu.:11.000   3rd Qu.: 0.1429  
##  Max.   :2014   Max.   :22.000   Max.   :17.0000  
##                                  NA's   :29
# First a plot of the eggs only from the breeding season:
par(mfrow = c(1,2), oma=c(0,0,2,0))

boxplot(eggs ~ year, data=ndat.egg, las = 2, ylab="Counts", ylim=c(-1,15), main="Mean counts per year")
title("Mean number of egg masses per year", outer=TRUE) 

boxplot(eggs ~ year, subset(ndat, eggs > 0), log = "y", ylab = "LOG(Counts)", las = 2, main="LOG(Mean counts per year when present)")

par(mfrow = c(1,2), oma=c(0,0,2,0))

boxplot(eggs ~ year, data=ndat, las = 2, ylab="Counts", ylim=c(-1,15), main="Mean counts per year")
title("Mean number of egg masses per year", outer=TRUE) 

boxplot(eggs ~ year, subset(ndat, eggs > 0), log = "y", ylab = "LOG(Counts)", las = 2, main="LOG(Mean counts per year when present)")

Code for Supplementary Figure S6

Variation in counts among lakes at Dusy Basin. Again, I am first showing LOG-transformed values! Also: I am splitting up again and showing egg data only for breeding season!

par(mfrow = c(3,1), mar = c(3,3,2,1), oma = c(2,2,0,0))

boxplot(larvae ~ site, subset(ndat, larvae > 0), log = "y", main = "Tadpoles (when present)", las = 2)
boxplot(subadults ~ site, subset(ndat, subadults > 0), log = "y", main = "Subadults (when present)", las = 2)
boxplot(adults ~ site, subset(ndat, adults > 0), log = "y", main = "Adults (when present)", las = 1)
mtext("LOG(Mean counts per site)", 2, 0, TRUE)
mtext("Sites", 1, 0, TRUE)

# Now trying to plot this also for ALL lakes!

par(mfrow = c(1,3), mar = c(3,3,2,1), oma = c(2,2,0,0))
boxplot(larvae ~ site, data=ndat, main = "Tadpoles", las = 2)
boxplot(subadults ~ site, data=ndat, main = "Subadults", las = 2)
boxplot(adults ~ site, data=ndat, main = "Adults", las = 1)
mtext("Mean counts per site", 2, 0, TRUE)
mtext("Sites", 1, 0, TRUE)

# And now an individual plot for the eggs, again first only for the breeding season:
quartz("", 12, 6)
par(mfrow = c(1,2), mar = c(3,3,2,1), oma = c(2,2,0,0))
boxplot(eggs ~ site, data=ndat.egg, main = "Egg masses", las = 2)
boxplot(eggs ~ site, subset(ndat, eggs > 0), log = "y", main = "LOG(Egg masses)", las = 2)
mtext("Sites", 1, 0, TRUE, cex=1.2)
mtext("Counts", 2, 0, TRUE, cex=1.2)

# And then for all seasons:
quartz("", 12, 6)
par(mfrow = c(1,2), mar = c(3,3,2,1), oma = c(2,2,0,0))
boxplot(eggs ~ site, data=ndat, main = "Egg masses", las = 2)
boxplot(eggs ~ site, subset(ndat, eggs > 0), log = "y", main = "LOG(Egg masses when present)", las = 2)
mtext("Sites", 1, 0, TRUE, cex=1.2)
mtext("Counts", 2, 0, TRUE, cex=1.2)

Code for Supplementary Figure S8

How does fish presence affect frog abundance at different life stages?

# Lake breeding quality

eg_max_y <- tapply(data$eggs, list(data$year, data$site), max, na.rm = TRUE)
breed_eggs <- apply(eg_max_y, 2, mean, na.rm = TRUE) / sum(apply(eg_max_y, 2, mean, na.rm = TRUE))

la_max_y <- tapply(data$larvae, list(data$year, data$site), max, na.rm = TRUE)
breed_larv <- apply(la_max_y, 2, mean, na.rm = TRUE) / sum(apply(la_max_y, 2, mean, na.rm = TRUE))

#quartz("", 12, 5)
par(mfrow = c(1,2), mar = c(5,4,1,1))
barplot(rbind(breed_eggs, breed_larv), beside = TRUE, las = 2)
mtext("Average number of egg masses / larvae (standardized)", 2, 3, las = 0)
plot(breed_larv[-c(9, 11, 13)] ~ apply(min_surf, 2, mean, na.rm = TRUE), pch = 16, log = "x", xlab = "Lake surface", xlim = c(0.3, 150000))
points(breed_eggs[-c(9,11,13)] ~ apply(min_surf, 2, mean, na.rm = TRUE), pch = 17, log = "x")
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "log" is not a
## graphical parameter
text(apply(min_surf, 2, mean, na.rm = TRUE), breed_larv[-c(9,11,13)], paste("Lake n°", names(breed_eggs[-c(9,11,13)])), pos = c(3, 4, 4, 4, 4, 4,3 , 4,3 ,4 , 4))

# Hypothesis 1: Best lakes are larger lakes / but fish eat frogs

fish <- table(data$fish, data$site, data$year)
test <- fish[2,,] / (fish[1,,] + fish[2,,])
test <- apply(test, 1, mean, na.rm = TRUE)

fish2 <- table(data$fish, data$site)
test2 <- fish2[2,] / (fish2[1,] + fish2[2,])

par(mfrow = c(1,2), mar = c(5,4,1,1), las = 1)
boxplot(min_surf+1, log = "y", las = 2, ylab = "Minimal lake surface")
#barplot(t(test), beside = TRUE, ylab = "Proportion of observed fish")
barplot(rbind(test, test2), beside = TRUE, ylab = "Proportion of observed fish", las =2)

lake_quality <- data.frame(breed_larv, breed_eggs, lake = as.numeric(names(breed_larv)))
prop_fish <- data.frame(prop_fish = test, lake = as.numeric(names(test)))
surface <- data.frame(surf.min = apply(min_surf, 2, mean, na.rm = TRUE), surf.mean = apply(mean_surf, 2, mean, na.rm = TRUE), surf.max = apply(max_surf, 2, max, na.rm = TRUE), lake = as.numeric(colnames(mean_surf)))
lake_quality <- merge(lake_quality, prop_fish)
lake_quality <- merge(lake_quality, surface)
lake_quality <- lake_quality[order(lake_quality$lake),]
lake_quality$fish <- "fishless"
    lake_quality$fish[lake_quality$lake %in% c(1,3)] <- "fishes"

# Lake 1, lake 3, early (lake 20 - fish died in 2000)

par(las = 1, plt = c(0,0.7, 0, 1), oma = c(4,5,1,1))
plot(breed_larv ~ surf.max, subset(lake_quality, fish == "fishless"), pch = 16, ylab = "", xlab = "", ylim = c(0, 0.5), xlim = c(0, 2800))
points(breed_eggs ~ surf.max, lake_quality, pch = 17)
res <- lm(breed_larv ~ surf.max + I(surf.max^2), subset(lake_quality, fish == "fishless"))
res <- lm(breed_larv ~ I(surf.max^2), subset(lake_quality, fish == "fishless"))
pred <- predict(res, newdata = data.frame(surf.max = c(0:3000)))
points(pred ~ c(0:3000), type = "l")
text(lake_quality$surf.max, lake_quality$breed_larv, lake_quality$lake, pos = c(rep(3, 7), 4, rep(3, 3)))
points(breed_larv ~ surf.max, subset(lake_quality, fish == "fishes"), pch = 1)
par(las = 1, plt = c(0.72,1, 0, 1), new = TRUE)
plot(breed_larv ~ surf.max, lake_quality, pch = 1, ylab = "", xlab = "", ylim = c(0, 0.5), xlim = c(50500, 51500), yaxt = "n", xaxt = "n")
axis(1, at = c(50500, 51000, 51500))
text(lake_quality$surf.max, lake_quality$breed_larv, lake_quality$lake, pos = 3)
mtext("Average number of tadpoles (standardized)", 2, 3, TRUE, las = 0)
mtext(expression(paste("Lake surface [", m^2, "]", sep = "")), 1, 3, TRUE)

# And as Steve Beissinger suggested let's try to also plot number of tadpoles in logarithm:

par(las = 1, plt = c(0,0.7, 0, 1), oma = c(4,5,1,1))

breed_larv <- log(breed_larv)
plot(breed_larv ~ surf.max, subset(lake_quality, fish == "fishless"), pch = 16, ylab = "", xlab = "", ylim = c(0, 0.5), xlim = c(0, 2800))
# points(breed_eggs ~ surf.max, lake_quality, pch = 17)
res <- lm(breed_larv ~ surf.max + I(surf.max^2), subset(lake_quality, fish == "fishless"))
res <- lm(breed_larv ~ I(surf.max^2), subset(lake_quality, fish == "fishless"))
pred <- predict(res, newdata = data.frame(surf.max = c(0:3000)))
points(pred ~ c(0:3000), type = "l")
text(lake_quality$surf.max, lake_quality$breed_larv, lake_quality$lake, pos = c(rep(3, 7), 4, rep(3, 3)))
points(breed_larv ~ surf.max, subset(lake_quality, fish == "fishes"), pch = 1)
par(las = 1, plt = c(0.72,1, 0, 1), new = TRUE)
plot(breed_larv ~ surf.max, lake_quality, pch = 1, ylab = "", xlab = "", ylim = c(0, 0.5), xlim = c(50500, 51500), yaxt = "n", xaxt = "n")
axis(1, at = c(50500, 51000, 51500))
text(lake_quality$surf.max, lake_quality$breed_larv, lake_quality$lake, pos = 3)
mtext("Average number of tadpoles (standardized)", 2, 3, TRUE, las = 0)
mtext(expression(paste("Lake surface [", m^2, "]", sep = "")), 1, 3, TRUE)

From here on, I am not running this code for publication anymore. It was used for data exploration. However, if you want to see all the analyses we did, turn eval=TRUE and it will run like a charm.

Now we are moving on to environmental factors!

Hypothesis: The years with low snow cover, or drought, are when lakes have dried and when tadpoles have died.  

Investigate for all different life stages:

wc <- tapply(dsnow$watcont, dsnow$year3, max, na.rm = TRUE)
wc <- data.frame(year = as.numeric(rownames(wc)), wc.max = wc, wc.mean = tapply(dsnow$watcont, dsnow$year3, mean, na.rm = TRUE))


# Season length estimated from figure with daily data
# ndat = new data, max per site and year

no_data2004 <- ndat[1:14,]
    no_data2004$site <- unique(ndat$site)
    no_data2004$year <- 2004
    no_data2004[,3:9] <- NA
no_data1996 <- no_data2004; no_data1996$year <- 1996
no_data1995 <- no_data2004; no_data1995$year <- 1995
no_data1994 <- no_data2004; no_data1994$year <- 1994
no_data1993 <- no_data2004; no_data1993$year <- 1993
no_data1992 <- no_data2004; no_data1992$year <- 1992

ndat <- rbind(ndat, no_data1992, no_data1993, no_data1994, no_data1995, no_data1996, no_data2004)

data2 <- merge(ndat, sle, all = TRUE)
data2 <- merge(data2, wc, all = TRUE)


# What environmental factor explains eggs abundance

res <- step(lmer(eggs ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res
## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                  13 -341.78 709.57                         
## (1 | site)          0   12 -359.35 742.71 35.145  1  3.061e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                     Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## sum.years:wc.max             1   0.133   0.133     1 128.42  0.0208
## sum.years:win.years          2   2.253   2.253     1 129.41  0.3567
## win.years:wc.max             3   0.163   0.163     1 130.69  0.0260
## win.years:winter             4  18.239  18.239     1 131.01  2.9243
## sum.years:winter             5   6.844   6.844     1 132.34  1.0824
## winter:wc.max                6  21.473  21.473     1 133.53  3.3928
## wc.max                       7  12.900  12.900     1 134.91  1.9992
## winter                       8   9.120   9.120     1 135.07  1.4002
## sum.years                    0 115.062 115.062     1 136.02 17.6067
## win.years                    0  53.193  53.193     1 136.04  8.1396
##                        Pr(>F)    
## sum.years:wc.max     0.885482    
## sum.years:win.years  0.551392    
## win.years:wc.max     0.872230    
## win.years:winter     0.089622 .  
## sum.years:winter     0.300056    
## winter:wc.max        0.067702 .  
## wc.max               0.159687    
## winter               0.238764    
## sum.years           4.873e-05 ***
## win.years            0.005008 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## eggs ~ sum.years + win.years + (1 | site)
res <- lmer(eggs ~ sum.years + win.years + ( 1|site), subset(data2, year < 2010)); summary(res); anova(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: eggs ~ sum.years + win.years + (1 | site)
##    Data: subset(data2, year < 2010)
## 
## REML criterion at convergence: 731
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3686 -0.2727 -0.0549  0.1046  6.8095 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 3.340    1.828   
##  Residual             6.535    2.556   
## Number of obs: 152, groups:  site, 14
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   14.465      3.531 141.600   4.096 7.03e-05 ***
## sum.years    -15.630      3.725 136.021  -4.196 4.87e-05 ***
## win.years    -11.092      3.888 136.036  -2.853  0.00501 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) sm.yrs
## sum.years -0.904       
## win.years -0.911  0.685
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## sum.years 115.062 115.062     1 136.02 17.6067 4.873e-05 ***
## win.years  53.193  53.193     1 136.04  8.1396  0.005008 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
a <- 100
pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), win.years = rep(seq(0.2, 0.8, length.out = a), times = a), site = rep(2, a*a)))

par(mar = c(4,4,1,1), las = 1)
image(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), xlab = "Summer length", ylab = "Winter length")
contour(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)

# What environmental factor explains larvae abundance

res <- step(lmer(larvae ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res
## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                  13 -963.25 1952.5                         
## (1 | site)          0   12 -992.64 2009.3 58.768  1  1.774e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                     Eliminated Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)
## sum.years:winter             1     12      12     1 128.38  0.0003 0.98624
## win.years:winter             2   1346    1346     1 129.10  0.0329 0.85645
## winter:wc.max                3   9065    9065     1 130.37  0.2229 0.63761
## win.years:wc.max             4  13515   13515     1 131.36  0.3343 0.56410
## sum.years:win.years          5   5733    5733     1 132.71  0.1425 0.70640
## winter                       6 144879  144879     1 133.13  3.6244 0.05910
## win.years                    0 224560  224560     1 134.21  5.5112 0.02036
## sum.years:wc.max             0 215705  215705     1 134.32  5.2939 0.02294
##                      
## sum.years:winter     
## win.years:winter     
## winter:wc.max        
## win.years:wc.max     
## sum.years:win.years  
## winter              .
## win.years           *
## sum.years:wc.max    *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## larvae ~ sum.years + win.years + wc.max + (1 | site) + sum.years:wc.max
res <- lmer(larvae ~ sum.years + win.years + ( 1|site), subset(data2, year < 2010)); summary(res); anova(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: larvae ~ sum.years + win.years + (1 | site)
##    Data: subset(data2, year < 2010)
## 
## REML criterion at convergence: 2042.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7947 -0.2347 -0.1101  0.1222  7.3332 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 34916    186.9   
##  Residual             41878    204.6   
## Number of obs: 152, groups:  site, 14
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)  
## (Intercept)    670.7      284.4  144.0   2.359   0.0197 *
## sum.years     -445.2      298.2  136.1  -1.493   0.1377  
## win.years     -638.5      311.3  136.1  -2.051   0.0422 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) sm.yrs
## sum.years -0.898       
## win.years -0.906  0.685
## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## sum.years  93354   93354     1 136.13  2.2292 0.13774  
## win.years 176226  176226     1 136.14  4.2080 0.04215 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), win.years = rep(seq(0.2, 0.8, length.out = a), times = a), site = rep(2, a*a)))
#quartz()
par(mar = c(4,4,1,1), las = 1)
image(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), xlab = "Summer length", ylab = "Winter length")
contour(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)

res <- lmer(larvae ~ sum.years*wc.max + ( 1|site), subset(data2, year < 2010)); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: larvae ~ sum.years * wc.max + (1 | site)
##    Data: subset(data2, year < 2010)
## 
## REML criterion at convergence: 2044.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8225 -0.2484 -0.1318  0.1175  7.1305 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 36429    190.9   
##  Residual             42066    205.1   
## Number of obs: 152, groups:  site, 14
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)        752.40     500.29   138.29   1.504   0.1349  
## sum.years        -1355.69     919.61   135.36  -1.474   0.1427  
## wc.max             -18.30      11.56   135.37  -1.583   0.1158  
## sum.years:wc.max    41.71      22.80   135.34   1.829   0.0696 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sm.yrs wc.max
## sum.years   -0.983              
## wc.max      -0.945  0.969       
## sm.yrs:wc.m  0.869 -0.927 -0.973
pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), wc.max = rep(seq(10, 60, length.out = a), times = a), site = rep(2, a*a)))
image(matrix(pred, nrow = a), xlab = "Summer length", ylab = "Water max")
contour(matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)

# What environmental factor explains juvenile abundance


res <- step(lmer(juveniles ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res
## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                  13 -756.14 1538.3                         
## (1 | site)          0   12 -777.45 1578.9 42.611  1   6.68e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                     Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## win.years:winter             1     53      53     1 128.26  0.0239
## sum.years:wc.max             2   1496    1496     1 128.96  0.6744
## winter:wc.max                3   2557    2557     1 130.11  1.1547
## sum.years:win.years          0  44577   44577     1 131.89 20.1195
## sum.years:winter             0  18088   18088     1 131.35  8.1639
## win.years:wc.max             0  31898   31898     1 131.42 14.3967
##                        Pr(>F)    
## win.years:winter     0.877434    
## sum.years:wc.max     0.413051    
## winter:wc.max        0.284552    
## sum.years:win.years 1.566e-05 ***
## sum.years:winter     0.004970 ** 
## win.years:wc.max     0.000225 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## juveniles ~ sum.years + win.years + winter + wc.max + (1 | site) + 
##     sum.years:win.years + sum.years:winter + win.years:wc.max
    res <- lmer(juveniles ~ sum.years * win.years + ( 1|site), subset(data2, year < 2010)); summary(res); anova(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: juveniles ~ sum.years * win.years + (1 | site)
##    Data: subset(data2, year < 2010)
## 
## REML criterion at convergence: 1604
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6173 -0.2714 -0.1482  0.1142  5.8718 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 1575     39.69   
##  Residual             2513     50.13   
## Number of obs: 152, groups:  site, 14
## 
## Fixed effects:
##                     Estimate Std. Error     df t value Pr(>|t|)
## (Intercept)            320.0      234.0  136.5   1.368    0.174
## sum.years             -470.4      424.9  136.0  -1.107    0.270
## win.years             -563.6      419.8  136.0  -1.342    0.182
## sum.years:win.years    897.7      780.5  136.0   1.150    0.252
## 
## Correlation of Fixed Effects:
##             (Intr) sm.yrs wn.yrs
## sum.years   -0.987              
## win.years   -0.988  0.990       
## sm.yrs:wn.y  0.955 -0.985 -0.983
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## sum.years           3080.3  3080.3     1 135.98  1.2255 0.2702
## win.years           4528.8  4528.8     1 136.03  1.8018 0.1817
## sum.years:win.years 3325.4  3325.4     1 136.05  1.3230 0.2521
pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), win.years = rep(seq(0.2, 0.8, length.out = a), times = a), site = rep(2, a*a)))

par(mar = c(4,4,1,1), las = 1)
image(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), xlab = "Summer length", ylab = "Winter length")
contour(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)

# What environmental factor explains subadults abundance

res <- step(lmer(subadults ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res
## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                  13 -578.76 1183.5                         
## (1 | site)          0   12 -588.81 1201.6 20.104  1  7.336e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                     Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## sum.years:win.years          0 1418.48 1418.48     1 129.00  7.4538
## sum.years:winter             0  896.33  896.33     1 128.65  4.7100
## sum.years:wc.max             0  788.27  788.27     1 128.60  4.1422
## win.years:winter             0  874.24  874.24     1 128.64  4.5939
## win.years:wc.max             0 1256.32 1256.32     1 128.76  6.6017
## winter:wc.max                0 1704.74 1704.74     1 128.12  8.9580
##                       Pr(>F)   
## sum.years:win.years 0.007216 **
## sum.years:winter    0.031824 * 
## sum.years:wc.max    0.043882 * 
## win.years:winter    0.033967 * 
## win.years:wc.max    0.011328 * 
## winter:wc.max       0.003315 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## subadults ~ (sum.years + win.years + winter + wc.max)^2 + (1 | 
##     site)
    res <- lmer(subadults ~ sum.years * win.years + ( 1|site), subset(data2, year < 2010)); summary(res); anova(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: subadults ~ sum.years * win.years + (1 | site)
##    Data: subset(data2, year < 2010)
## 
## REML criterion at convergence: 1231.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4208 -0.5041 -0.0731  0.0982  6.3068 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept)  69.41    8.331  
##  Residual             212.85   14.589  
## Number of obs: 152, groups:  site, 14
## 
## Fixed effects:
##                     Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)          -102.88      67.98  136.87  -1.513   0.1325  
## sum.years             219.77     123.47  136.63   1.780   0.0773 .
## win.years             180.60     122.00  136.71   1.480   0.1411  
## sum.years:win.years  -362.92     226.80  136.73  -1.600   0.1119  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sm.yrs wn.yrs
## sum.years   -0.987              
## win.years   -0.989  0.990       
## sm.yrs:wn.y  0.955 -0.985 -0.983
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## sum.years           674.32  674.32     1 136.63  3.1680 0.07732 .
## win.years           466.44  466.44     1 136.71  2.1914 0.14109  
## sum.years:win.years 545.04  545.04     1 136.73  2.5607 0.11186  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), win.years = rep(seq(0.2, 0.8, length.out = a), times = a), site = rep(2, a*a)))

par(mar = c(4,4,1,1), las = 1)
image(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), xlab = "Summer length", ylab = "Winter length")
contour(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)

# What environmental factor explains adults abundance

res <- step(lmer(adults ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res
## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                  13 -544.36 1114.7                         
## (1 | site)          0   12 -588.25 1200.5 87.781  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                     Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## winter:wc.max                1  102.48  102.48     1 128.02  0.9872
## sum.years:winter             2   36.16   36.16     1 129.21  0.3485
## sum.years:wc.max             3   14.40   14.40     1 130.12  0.1395
## sum.years:win.years          4   68.68   68.68     1 131.33  0.6698
## sum.years                    5   15.99   15.99     1 132.24  0.1564
## win.years:winter             0  473.20  473.20     1 133.02  4.6599
## win.years:wc.max             0 1830.34 1830.34     1 133.23 18.0244
##                        Pr(>F)    
## winter:wc.max         0.32230    
## sum.years:winter      0.55601    
## sum.years:wc.max      0.70936    
## sum.years:win.years   0.41459    
## sum.years             0.69315    
## win.years:winter      0.03267 *  
## win.years:wc.max    4.062e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## adults ~ win.years + winter + wc.max + (1 | site) + win.years:winter + 
##     win.years:wc.max
# What environmental factor explains adults and subadults abundance

res <- step(lmer(frogs ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res
## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                  13 -640.94 1307.9                         
## (1 | site)          0   12 -671.15 1366.3 60.424  1  7.648e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                     Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## sum.years:wc.max             1 1555.6  1555.6     1 128.20  3.6733
## sum.years:winter             2  325.6   325.6     1 128.97  0.7545
## win.years:winter             3  270.1   270.1     1 129.97  0.6274
## sum.years:win.years          0 1890.9  1890.9     1 131.34  4.4069
## win.years:wc.max             0 5852.7  5852.7     1 131.09 13.6404
## winter:wc.max                0 1873.0  1873.0     1 131.37  4.3653
##                        Pr(>F)    
## sum.years:wc.max    0.0575177 .  
## sum.years:winter    0.3866705    
## win.years:winter    0.4297400    
## sum.years:win.years 0.0377068 *  
## win.years:wc.max    0.0003238 ***
## winter:wc.max       0.0386071 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## frogs ~ sum.years + win.years + winter + wc.max + (1 | site) + 
##     sum.years:win.years + win.years:wc.max + winter:wc.max

Drought and season length investigations

How season length affects drying out of pools and lakes:

table(dry$dryness, dry$lake)
##      
##        2  5  6  7  8  9 10 11
##   dry  8  0  8  0 11  7 10  7
##   low  2  3  0  5  0  1  0  0
##   wet  4 11  6  9  3  6  4  2
test <- merge(sle, dry)

res <- lmer(winter ~ dryness + (1|lake), test); anova(res)
## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq  Mean Sq NumDF DenDF F value   Pr(>F)   
## dryness 0.070291 0.035146     2   104  4.8678 0.009531 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res <- lmer(sum.years ~ dryness + (1|lake), test); anova(res)
## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## dryness 0.05313 0.026565     2   104  4.2857 0.01627 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(1,2), las = 1, mar = c(3,4,1,1), oma = c(0,0,2,1))
boxplot(winter ~ dryness, subset(test, lake == 2), ylab = "Winter length")
    res <- lm(winter ~ dryness, subset(test, lake == 2)); anova(res)
## Analysis of Variance Table
## 
## Response: winter
##           Df   Sum Sq   Mean Sq F value  Pr(>F)  
## dryness    2 0.050474 0.0252369  4.9102 0.02992 *
## Residuals 11 0.056537 0.0051397                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    res2 <- lmer(winter ~ dryness + (1|lake), test); anova(res2)
## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq  Mean Sq NumDF DenDF F value   Pr(>F)   
## dryness 0.070291 0.035146     2   104  4.8678 0.009531 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    legend("topleft", c(paste(p.trans(anova(res)$Pr[1]), "(lake 2)"), paste(p.trans(anova(res2)$Pr[1]), "(all lakes)")), bty = "n")
boxplot(sum.years ~ dryness, subset(test, lake == 2), ylab = "Summer length")
    res <- lm(sum.years ~ dryness, subset(test, lake == 2)); anova(res)
## Analysis of Variance Table
## 
## Response: sum.years
##           Df   Sum Sq   Mean Sq F value Pr(>F)
## dryness    2 0.024695 0.0123477  2.0327 0.1773
## Residuals 11 0.066821 0.0060746
    res2 <- lmer(sum.years ~ dryness + (1|lake), test); anova(res2)
## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## dryness 0.05313 0.026565     2   104  4.2857 0.01627 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    legend("bottomleft", c(paste(p.trans(anova(res)$Pr[1]), "(lake 2)"), paste(p.trans(anova(res2)$Pr[1]), "(all lakes)")), bty = "n")
mtext("Lake 2 example", 3, 0, font = 2, out = TRUE)

Relationship between drought (in drying lakes) and snow pack and season lengths.

# Combined files year of drought versus year after drought

test <- merge(sle, dry)
drought <- merge(dry, data)
drought <- subset(drought, season == "b")
dry2 <- test
dry2$year <- I(dry2$year + 1)
drought2 <- merge(dry2, data)
drought2 <- subset(drought2, season == "b")

dry$site <- dry$lake
drought <- merge(dry, data2)


dry2 <- dry
dry2$year <- I(dry$year + 1)
drought2 <- merge(dry2, data2)


# Data structure

table(drought$year, drought$site)
##       
##        2 5 6 7 8 9 10 11
##   1997 1 1 1 1 1 1  1  1
##   1998 1 1 1 1 1 1  1  1
##   1999 1 1 1 1 1 1  1  1
##   2000 1 1 1 1 1 1  1  1
##   2001 1 1 1 1 1 1  1  1
##   2002 1 1 1 1 1 1  1  1
##   2003 1 1 1 1 1 1  1  1
##   2004 1 1 1 1 1 1  1  1
##   2005 1 1 1 1 1 1  1  1
##   2006 1 1 1 1 1 1  1  1
##   2007 1 1 1 1 1 1  1  1
##   2008 1 1 1 1 1 1  1  1
##   2009 1 1 1 1 1 1  1  1
##   2010 1 1 1 1 1 1  1  1
table(drought$site, drought$dryness)
##     
##      dry low wet
##   2    8   2   4
##   5    0   3  11
##   6    8   0   6
##   7    0   5   9
##   8   11   0   3
##   9    7   1   6
##   10  10   0   4
##   11   7   0   2
# Combined files year of drought versus year after drought

par(mfrow = c(2,4), mar = c(1,2,2,1), oma = c(4,3,1,0))
boxplot(I(eggs+1) ~ dryness, drought, log = "y", main = "Egg masses", las = 1)
    res <- lmer(log(eggs+1) ~ dryness + (1|site), drought); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(eggs + 1) ~ dryness + (1 | site)
##    Data: drought
## 
## REML criterion at convergence: 175.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.95636 -0.27087 -0.06274 -0.02951  2.80811 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.1438   0.3793  
##  Residual             0.3175   0.5634  
## Number of obs: 94, groups:  site, 8
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)  0.28249    0.16258  9.76009   1.738    0.114
## drynesslow   0.07577    0.25000 90.98213   0.303    0.763
## drynesswet   0.11480    0.14308 90.99677   0.802    0.424
## 
## Correlation of Fixed Effects:
##            (Intr) drynssl
## drynesslow -0.284        
## drynesswet -0.421  0.433
    legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Pop. size (during drought)", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ dryness, drought, log = "y", main = "Larvae", las = 1)
    res <- lmer(log(larvae +1) ~ dryness + (1|site), drought); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(larvae + 1) ~ dryness + (1 | site)
##    Data: drought
## 
## REML criterion at convergence: 382
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.03642 -0.57303 -0.05086  0.57345  2.43186 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 3.666    1.915   
##  Residual             2.851    1.688   
## Number of obs: 94, groups:  site, 8
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)   2.1603     0.7325  8.0076   2.949   0.0184 *
## drynesslow    0.7316     0.7663 88.2544   0.955   0.3423  
## drynesswet    0.1942     0.4393 88.4867   0.442   0.6595  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) drynssl
## drynesslow -0.198        
## drynesswet -0.288  0.456
    legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(juveniles+1) ~ dryness, drought, log = "y", main = "Juveniles", las = 1)
    res <- lmer(log(juveniles +1) ~ dryness + (1|site), drought); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(juveniles + 1) ~ dryness + (1 | site)
##    Data: drought
## 
## REML criterion at convergence: 305
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6818 -0.3661 -0.1764  0.1905  2.7698 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 1.016    1.008   
##  Residual             1.264    1.124   
## Number of obs: 94, groups:  site, 8
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)  0.72073    0.40165  8.50316   1.794    0.108
## drynesslow   0.27373    0.50614 89.75381   0.541    0.590
## drynesswet   0.04649    0.28998 89.99626   0.160    0.873
## 
## Correlation of Fixed Effects:
##            (Intr) drynssl
## drynesslow -0.236        
## drynesswet -0.346  0.448
    legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(frogs+1) ~ dryness, drought, log = "y", main = "Adults and subadults", las = 1)
    res <- lmer(log(frogs +1) ~ dryness + (1|site), drought); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(frogs + 1) ~ dryness + (1 | site)
##    Data: drought
## 
## REML criterion at convergence: 267.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.10960 -0.51168 -0.02874  0.57383  2.39248 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.5294   0.7276  
##  Residual             0.8540   0.9241  
## Number of obs: 94, groups:  site, 8
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)   1.3780     0.2986  9.0854   4.614  0.00123 **
## drynesslow    0.4244     0.4136 90.4971   1.026  0.30760   
## drynesswet    0.4616     0.2368 90.6815   1.949  0.05438 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) drynssl
## drynesslow -0.258        
## drynesswet -0.380  0.442
    legend("topleft", p.trans(anova(res)$P), bty = "n")

boxplot(I(eggs+1) ~ dryness, drought2, log = "y", main = "", las = 1)
    res <- lmer(log(eggs+1) ~ dryness + (1|site), drought2); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(eggs + 1) ~ dryness + (1 | site)
##    Data: drought2
## 
## REML criterion at convergence: 172.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.11230 -0.42444  0.01234  0.11625  2.83040 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.1360   0.3688  
##  Residual             0.3115   0.5581  
## Number of obs: 93, groups:  site, 8
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)   0.2013     0.1602 10.5685   1.256   0.2361  
## drynesslow    0.2518     0.2195 89.9930   1.147   0.2543  
## drynesswet    0.2625     0.1399 89.7020   1.876   0.0639 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) drynssl
## drynesslow -0.317        
## drynesswet -0.430  0.453
    legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Pop. size (after drought)", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ dryness, drought2, log = "y", main = "", las = 1)
    res <- lmer(log(larvae +1) ~ dryness + (1|site), drought2); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(larvae + 1) ~ dryness + (1 | site)
##    Data: drought2
## 
## REML criterion at convergence: 386.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.95250 -0.59904 -0.04415  0.63750  2.32732 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 3.246    1.802   
##  Residual             3.166    1.779   
## Number of obs: 93, groups:  site, 8
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)   2.0328     0.7043  8.5488   2.886    0.019 *
## drynesslow    1.1589     0.7140 87.9479   1.623    0.108  
## drynesswet    0.2407     0.4535 87.2611   0.531    0.597  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) drynssl
## drynesslow -0.238        
## drynesswet -0.318  0.471
    legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(juveniles+1) ~ dryness, drought2, log = "y", main = "", las = 1)
    res <- lmer(log(juveniles +1) ~ dryness + (1|site), drought2); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(juveniles + 1) ~ dryness + (1 | site)
##    Data: drought2
## 
## REML criterion at convergence: 300.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6295 -0.3508 -0.1245  0.1737  2.8947 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.881    0.9386  
##  Residual             1.253    1.1196  
## Number of obs: 93, groups:  site, 8
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)   0.5977     0.3815  9.0381   1.567    0.151
## drynesslow    0.6498     0.4459 89.1122   1.457    0.149
## drynesswet    0.2507     0.2836 88.4128   0.884    0.379
## 
## Correlation of Fixed Effects:
##            (Intr) drynssl
## drynesslow -0.273        
## drynesswet -0.367  0.464
    legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(frogs+1) ~ dryness, drought2, log = "y", main = "", las = 1)
    res <- lmer(log(frogs +1) ~ dryness + (1|site), drought2); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(frogs + 1) ~ dryness + (1 | site)
##    Data: drought2
## 
## REML criterion at convergence: 272.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.88822 -0.43647 -0.03915  0.60032  2.47827 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.4261   0.6528  
##  Residual             0.9492   0.9743  
## Number of obs: 93, groups:  site, 8
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)   1.1324     0.2824 10.2815   4.011  0.00234 **
## drynesslow    0.7393     0.3835 89.9802   1.928  0.05705 . 
## drynesswet    0.5953     0.2444 89.6422   2.436  0.01684 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) drynssl
## drynesslow -0.315        
## drynesswet -0.426  0.454
    legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Dry versus not dry lake", side = 1, line = 2, outer = TRUE)

Same figure for lake 2 only

par(mfrow = c(2,4), mar = c(1,2,2,1), oma = c(4,3,1,0))
boxplot(I(eggs+1) ~ dryness, subset(drought, site == 2), log = "y", main = "Egg masses", las = 1)
    res <- lm(log(eggs+1) ~ dryness, subset(drought, site == 2)); summary(res)
## 
## Call:
## lm(formula = log(eggs + 1) ~ dryness, data = subset(drought, 
##     site == 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2184 -1.2043 -0.2740  0.9563  1.3465 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  1.21844    0.44934   2.712   0.0219 *
## drynesslow  -0.38932    0.95320  -0.408   0.6916  
## drynesswet  -0.01419    0.74515  -0.019   0.9852  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.189 on 10 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.0174, Adjusted R-squared:  -0.1791 
## F-statistic: 0.08854 on 2 and 10 DF,  p-value: 0.916
    legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
mtext("Pop. size (during subset(drought, site == 2))", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ dryness, subset(drought, site == 2), log = "y", main = "Larvae", las = 1)
    res <- lm(log(larvae +1) ~ dryness, subset(drought, site == 2)); summary(res)
## 
## Call:
## lm(formula = log(larvae + 1) ~ dryness, data = subset(drought, 
##     site == 2))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.635 -2.820  1.220  1.964  3.016 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    4.635      1.146   4.044  0.00235 **
## drynesslow    -1.619      2.431  -0.666  0.52056   
## drynesswet     1.026      1.901   0.540  0.60127   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.032 on 10 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.09267,    Adjusted R-squared:  -0.08879 
## F-statistic: 0.5107 on 2 and 10 DF,  p-value: 0.6149
    legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(juveniles+1) ~ dryness, subset(drought, site == 2), log = "y", main = "Juveniles", las = 1)
    res <- lm(log(juveniles +1) ~ dryness, subset(drought, site == 2)); summary(res)
## 
## Call:
## lm(formula = log(juveniles + 1) ~ dryness, data = subset(drought, 
##     site == 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0636 -0.6287 -0.3102  0.8974  2.0636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.6287     0.4526   1.389    0.195
## drynesslow    1.4349     0.9602   1.494    0.166
## drynesswet   -0.1685     0.7506  -0.225    0.827
## 
## Residual standard error: 1.198 on 10 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2121, Adjusted R-squared:  0.05453 
## F-statistic: 1.346 on 2 and 10 DF,  p-value: 0.3036
    legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(frogs+1) ~ dryness, subset(drought, site == 2), log = "y", main = "Adults and subadults", las = 1)
    res <- lm(log(frogs +1) ~ dryness, subset(drought, site == 2)); summary(res)
## 
## Call:
## lm(formula = log(frogs + 1) ~ dryness, data = subset(drought, 
##     site == 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5423 -0.6908  0.4268  0.8693  2.5423 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  2.48253    0.58094   4.273  0.00163 **
## drynesslow   0.05973    1.23237   0.048  0.96230   
## drynesswet   0.96780    0.96339   1.005  0.33879   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.537 on 10 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.09669,    Adjusted R-squared:  -0.08398 
## F-statistic: 0.5352 on 2 and 10 DF,  p-value: 0.6014
    legend("topleft", p.trans(anova(res)$P[1]), bty = "n")

boxplot(I(eggs+1) ~ dryness, subset(drought2, site == 2), log = "y", main = "", las = 1)
    res <- lm(log(eggs+1) ~ dryness, subset(drought2, site == 2)); summary(res)
## 
## Call:
## lm(formula = log(eggs + 1) ~ dryness, data = subset(drought2, 
##     site == 2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.28247 -0.76013 -0.06698  0.42889  1.50855 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.7601     0.3591   2.117   0.0603 .
## drynesslow    0.5223     0.8029   0.651   0.5300  
## drynesswet    1.3593     0.6875   1.977   0.0762 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.016 on 10 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.283,  Adjusted R-squared:  0.1396 
## F-statistic: 1.973 on 2 and 10 DF,  p-value: 0.1895
    legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
mtext("Pop. size (after subset(drought, site == 2))", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ dryness, subset(drought2, site == 2), log = "y", main = "", las = 1)
    res <- lm(log(larvae +1) ~ dryness, subset(drought2, site == 2)); summary(res)
## 
## Call:
## lm(formula = log(larvae + 1) ~ dryness, data = subset(drought2, 
##     site == 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0253 -3.2673  0.8187  1.9769  3.5987 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   4.0253     1.1312   3.558   0.0052 **
## drynesslow   -0.7581     2.5295  -0.300   0.7705   
## drynesswet    2.4870     2.1661   1.148   0.2776   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.2 on 10 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1423, Adjusted R-squared:  -0.02922 
## F-statistic: 0.8296 on 2 and 10 DF,  p-value: 0.4641
    legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(juveniles+1) ~ dryness, subset(drought2, site == 2), log = "y", main = "", las = 1)
    res <- lm(log(juveniles +1) ~ dryness, subset(drought2, site == 2)); summary(res)
## 
## Call:
## lm(formula = log(juveniles + 1) ~ dryness, data = subset(drought2, 
##     site == 2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.07709 -0.51713 -0.06515  0.02704  2.05004 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.5171     0.3773   1.371   0.2005  
## drynesslow   -0.5171     0.8437  -0.613   0.5536  
## drynesswet    1.5600     0.7225   2.159   0.0562 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.067 on 10 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3743, Adjusted R-squared:  0.2492 
## F-statistic: 2.991 on 2 and 10 DF,  p-value: 0.09589
    legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(frogs+1) ~ dryness, subset(drought2, site == 2), log = "y", main = "", las = 1)
    res <- lm(log(frogs +1) ~ dryness, subset(drought2, site == 2)); summary(res)
## 
## Call:
## lm(formula = log(frogs + 1) ~ dryness, data = subset(drought2, 
##     site == 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1841 -0.6057  0.4550  1.1636  1.1763 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   2.1841     0.4826   4.526   0.0011 **
## drynesslow   -1.0205     1.0791  -0.946   0.3666   
## drynesswet    1.9939     0.9241   2.158   0.0563 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.365 on 10 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4075, Adjusted R-squared:  0.289 
## F-statistic: 3.439 on 2 and 10 DF,  p-value: 0.07302
    legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
mtext("Dry versus not dry lake", side = 1, line = 2, outer = TRUE)

Time structured analyses

How does recruitment in lake 2 affect the whole population?

Check how the number of previous life stages correlates with numbers of later life stages.

# Functions for drawing

plot.regr.FUN <- function(x,y, alpha = 0.05, pch = 16, yaxt = "s", lab = names(x)) {
    plot(y ~ x, las = 1, pch = pch, yaxt = yaxt)
    text(x, y, lab, pos = 3, col = "darkgrey")
    res <- lm(y ~ x)
    if (anova(res)$"Pr(>F)"[1] <= alpha) {abline(res); legend("topright", p.trans(anova(res)$"Pr(>F)"[1]), bty = "n")}
    if (anova(res)$"Pr(>F)"[1] > alpha) {abline(h = mean(y, na.rm = TRUE), lty = 2, col = "darkgrey");legend("topright", "p > 0.05 [NS]", bty = "n")}
    }

time.FUN <- function(x,y, time = 5, ytitle = "") {
    plot.regr.FUN(x,y)
    mtext(ytitle, 2, 3, F, cex = 0.8)
    n <- length(x)
    for (i in 1:(time-1)) plot.regr.FUN(y = y[-c(1:i)], x = x[-c((n-(i-1)):n)], yaxt = "n")
}

# Temporal autocorrelation in lake 2

data3 <- subset(data2, site == 2 & year <= 2009)
x1 <- data3$eggs; names(x1) <- 1992:2009
x2 <- data3$larvae; names(x2) <- 1992:2009
x3 <- data3$juveniles; names(x3) <- 1992:2009
x4 <- data3$frogs; names(x4) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
time.FUN(x2, data3$larvae, ytitle = "Larvae")
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
time.FUN(x3, data3$juveniles, ytitle = "Juveniles")
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
time.FUN(x4, data3$frogs, ytitle = "Adults and subadults")
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
mtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Temporal autocorrelation in lake 2", 3, 1, TRUE, font = 2)

# Relationships between metrics

pairs(data3[,3:9])

pairs(data2[,10:14])

# Previous life stages effects in lake 2 on the whole bassin

data3 <- subset(data2, site == 2 & year <= 2009)
x1 <- tapply(data2$eggs, data2$year, mean, na.rm = TRUE)[4:21]
x2 <- tapply(data2$larvae, data2$year, mean, na.rm = TRUE)[4:21]
x3 <- tapply(data2$juveniles, data2$year, mean, na.rm = TRUE)[4:21]
x4 <- tapply(data2$frogs, data2$year, mean, na.rm = TRUE)[4:21]

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$eggs, x1, ytitle = "Egg masses")
time.FUN(data3$eggs, x2, ytitle = "Larvae")
time.FUN(data3$eggs, x3, ytitle = "Juveniles")
time.FUN(data3$eggs, x4, ytitle = "Adults and subadults")
mtext("Egg masses in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of eggs masses from lake 2 on the whole basin", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$larvae, x1, ytitle = "Egg masses")
time.FUN(data3$larvae, x2, ytitle = "Larvae")
time.FUN(data3$larvae, x3, ytitle = "Juveniles")
time.FUN(data3$larvae, x4, ytitle = "Adults and subadults")
mtext("Larvae in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of larvae from lake 2 on the whole basin", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$juveniles, x1, ytitle = "Egg masses")
time.FUN(data3$juveniles, x2, ytitle = "Larvae")
time.FUN(data3$juveniles, x3, ytitle = "Juveniles")
time.FUN(data3$juveniles, x4, ytitle = "Adults and subadults")
mtext("Juveniles in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of juveniles from lake 2 on the whole basin", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$frogs, x1, ytitle = "Egg masses")
time.FUN(data3$frogs, x2, ytitle = "Larvae")
time.FUN(data3$frogs, x3, ytitle = "Juveniles")
time.FUN(data3$frogs, x4, ytitle = "Adults and subadults")
mtext("Adults and subadults in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of adults and subadults from lake 2 on the whole basin", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$juveniles[-8], x1[-8], ytitle = "Egg masses")
time.FUN(data3$juveniles[-8], x2[-8], ytitle = "Larvae")
time.FUN(data3$juveniles[-8], x3[-8], ytitle = "Juveniles")
time.FUN(data3$juveniles[-8], x4[-8], ytitle = "Adults and subadults")
mtext("Juveniles in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of juveniles from lake 2 on the whole basin (without 1999)", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$frogs[-8], x1[-8], ytitle = "Egg masses")
time.FUN(data3$frogs[-8], x2[-8], ytitle = "Larvae")
time.FUN(data3$frogs[-8], x3[-8], ytitle = "Juveniles")
time.FUN(data3$frogs[-8], x4[-8], ytitle = "Adults and subadults")
mtext("Adults and subadults in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of adults and subadults from lake 2 on the whole basin (without 1999)", 3, 1, TRUE, font = 2)

The same but only in lake 2:

data3 <- subset(data2, site == 2 & year <= 2009)

x1 <- data3$eggs; names(x1) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
time.FUN(x1, data3$larvae, ytitle = "Larvae")
time.FUN(x1, data3$juveniles, ytitle = "Juveniles")
time.FUN(x1, data3$frogs, ytitle = "Adults and subadults")
mtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Number of egg masses in previous years", 3, 1, TRUE, font = 2)

x1 <- data3$larvae; names(x1) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")
time.FUN(x1, data3$larvae, ytitle = "Larvae")
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
time.FUN(x1, data3$juveniles, ytitle = "Juveniles")
time.FUN(x1, data3$frogs, ytitle = "Adults and subadults")
mtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Number of larvae in previous years", 3, 1, TRUE, font = 2)

x1 <- data3$juveniles; names(x1) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")
time.FUN(x1, data3$larvae, ytitle = "Larvae")
time.FUN(x1, data3$juveniles, ytitle = "Juveniles")
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
time.FUN(x1, data3$frogs, ytitle = "Adults and subadults")
mtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Number of juveniles in previous years", 3, 1, TRUE, font = 2)

x1 <- data3$frogs; names(x1) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")
time.FUN(x1, data3$larvae, ytitle = "Larvae")
time.FUN(x1, data3$juveniles, ytitle = "Juveniles")
time.FUN(x1, data3$frogs, ytitle = "Adults and subadults")
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
mtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Number of adults and subadults in previous years", 3, 1, TRUE, font = 2)

Lots of plots to figure out whether environmental data explains how many frogs were counted each summer:

# Proportion of occupied lakes as functions of seasonal data

YGF2 <- cbind(matrix(NA, 4,4), YGF[,1:7], rep(NA, 4), YGF[,8:12])
    colnames(YGF2) <- 1993:2009

par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF2[,-c(1:4)], beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = terrain.colors(4), ylab = "Proportion of occupied lakes")

x <- data3$winter[-1]; names(x) <- colnames(YGF2) <- 1993:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, YGF2[1,], ytitle = "Egg masses")
time.FUN(x, YGF2[2,], ytitle = "Larvae")
time.FUN(x, YGF2[3,], ytitle = "Juveniles")
time.FUN(x, YGF2[4,], ytitle = "Adults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Proportion of occupied lakes", 3, 1, TRUE, font = 2)

x <- data3$wc.mean[-1]; names(x) <- colnames(YGF2) <- 1993:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, YGF2[1,], ytitle = "Egg masses")
time.FUN(x, YGF2[2,], ytitle = "Larvae")
time.FUN(x, YGF2[3,], ytitle = "Juveniles")
time.FUN(x, YGF2[4,], ytitle = "Adults")
mtext("Average water content", 1, 1, T, cex = 0.8)
mtext("Proportion of occupied lakes", 3, 1, TRUE, font = 2)

x <- data3$sum.years[-1]; names(x) <- colnames(YGF2) <- 1993:2009
#quartz("", 12, 7)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, YGF2[1,], ytitle = "Egg masses")
time.FUN(x, YGF2[2,], ytitle = "Larvae")
time.FUN(x, YGF2[3,], ytitle = "Juveniles")
time.FUN(x, YGF2[4,], ytitle = "Adults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Proportion of occupied lakes", 3, 1, TRUE, font = 2)

# Winter length, water content (wc) and summer length in lake 2

data3 <- subset(data2, site == 2 & year <= 2009)

x <- data3$winter; names(x) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 2", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 2", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 2", 3, 1, TRUE, font = 2)

# Winter length, water content (wc) and summer length in lake 5

data3 <- subset(data2, site == 5 & year <= 2009)

x <- data3$winter; names(x) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 5", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 5", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 5", 3, 1, TRUE, font = 2)

# Winter length, water content (wc) and summer length in lake 7

data3 <- subset(data2, site == 7 & year <= 2009)

x <- data3$winter; names(x) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 7", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 7", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 7", 3, 1, TRUE, font = 2)

# Winter length, water content (wc) and summer length in lake 20

data3 <- subset(data2, site == 20 & year <= 2009)

x <- data3$winter; names(x) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 20", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 20", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 20", 3, 1, TRUE, font = 2)

# Winter length, water content (wc) and summer length in lake 4

data3 <- subset(data2, site == 4 & year <= 2009)

x <- data3$winter; names(x) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 4", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 4", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 4", 3, 1, TRUE, font = 2)

# Winter length, water content (wc) and summer length in lake 1

data3 <- subset(data2, site == 1 & year <= 2009)

x <- data3$winter; names(x) <- 1992:2009

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 1", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 1", 3, 1, TRUE, font = 2)

par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 1", 3, 1, TRUE, font = 2)

Population size vs. lake dried out or not

table(water$Year)
## 
## 2002 2003 2005 2006 
##   49   29   41   30
drought.prob <- data.frame(min_surf = as.numeric(min_surf), year = as.numeric(rep(rownames(min_surf), times = dim(min_surf)[2])), site = as.numeric(rep(colnames(min_surf), each = dim(min_surf)[1])))

#--------------------------------------------------------
# Lake dried (0), or not (1)
#--------------------------------------------------------
drought.prob$drought <- drought.prob$min_surf
    drought.prob$drought[drought.prob$drought > 0] <- 1
drought.prob2 <- drought.prob
    drought.prob2$year <- drought.prob2$year+1

#--------------------------------------------------------
# Match population summaries (data2) with drought prob
# The year during and after the drought
#--------------------------------------------------------
drought.prob <- merge(data2, drought.prob)
drought.prob2 <- merge(data2, drought.prob2)

#--------------------------------------------------------
# Compare populations between lakes that dried versus not dried
#--------------------------------------------------------
table(drought.prob$drought, drought.prob$site, drought.prob$year)
## , ,  = 2002
## 
##    
##     1 2 3 4 5 6 7 8 10 20 22
##   0 0 1 0 0 0 1 0 1  1  0  0
##   1 1 0 1 1 1 0 1 0  0  1  0
## 
## , ,  = 2003
## 
##    
##     1 2 3 4 5 6 7 8 10 20 22
##   0 0 0 0 0 0 1 0 1  1  0  0
##   1 1 1 1 1 1 0 1 0  0  0  1
## 
## , ,  = 2005
## 
##    
##     1 2 3 4 5 6 7 8 10 20 22
##   0 0 0 0 0 0 0 0 1  1  0  0
##   1 1 1 1 1 1 1 1 0  0  1  1
## 
## , ,  = 2006
## 
##    
##     1 2 3 4 5 6 7 8 10 20 22
##   0 0 1 0 0 0 0 0 0  0  0  0
##   1 1 0 1 1 1 1 1 1  1  1  0
par(mfrow = c(2,4), mar = c(1,2,2,1), oma = c(4,3,1,0))
boxplot(I(eggs+1) ~ drought, drought.prob, log = "y", main = "Egg masses", las = 1)
    res <- lmer(log(eggs+1) ~ drought + (1|site), drought.prob); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(eggs + 1) ~ drought + (1 | site)
##    Data: drought.prob
## 
## REML criterion at convergence: 94
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4924 -0.3973 -0.2604 -0.1466  2.2097 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.1691   0.4113  
##  Residual             0.4563   0.6755  
## Number of obs: 41, groups:  site, 11
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)  0.40387    0.27979 24.47525   1.444    0.162
## drought      0.06342    0.29716 37.59510   0.213    0.832
## 
## Correlation of Fixed Effects:
##         (Intr)
## drought -0.811
    legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Pop. size (during drought)", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ drought, drought.prob, log = "y", main = "Larvae", las = 1)
    res <- lmer(log(larvae +1) ~ drought + (1|site), drought.prob); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(larvae + 1) ~ drought + (1 | site)
##    Data: drought.prob
## 
## REML criterion at convergence: 147.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87044 -0.44350 -0.08483  0.48452  2.55453 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 6.919    2.6303  
##  Residual             0.959    0.9793  
## Number of obs: 41, groups:  site, 11
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   4.3667     0.9001 14.5264   4.851 0.000231 ***
## drought      -1.1922     0.5132 31.1300  -2.323 0.026879 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## drought -0.440
    legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(juveniles+1) ~ drought, drought.prob, log = "y", main = "Juveniles", las = 1)
    res <- lmer(log(juveniles +1) ~ drought + (1|site), drought.prob); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(juveniles + 1) ~ drought + (1 | site)
##    Data: drought.prob
## 
## REML criterion at convergence: 142.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.60978 -0.28228  0.03419  0.55254  1.49806 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 2.574    1.604   
##  Residual             1.130    1.063   
## Number of obs: 41, groups:  site, 11
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)   1.4202     0.6587 20.6056   2.156   0.0431 *
## drought       0.8475     0.5373 34.7451   1.577   0.1238  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## drought -0.628
    legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(frogs+1) ~ drought, drought.prob, log = "y", main = "Adults and subadults", las = 1)
    res <- lmer(log(frogs +1) ~ drought + (1|site), drought.prob); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(frogs + 1) ~ drought + (1 | site)
##    Data: drought.prob
## 
## REML criterion at convergence: 88
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.56119 -0.62143  0.05131  0.64305  1.51461 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 1.2259   1.1072  
##  Residual             0.2258   0.4751  
## Number of obs: 41, groups:  site, 11
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   2.1267     0.3920 15.8250   5.426 5.82e-05 ***
## drought       0.4606     0.2475 31.7976   1.861    0.072 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## drought -0.487
    legend("topleft", p.trans(anova(res)$P), bty = "n")

boxplot(I(eggs+1) ~ drought, drought.prob2, log = "y", main = "", las = 1)
    res <- lmer(log(eggs+1) ~ drought + (1|site), drought.prob2); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(eggs + 1) ~ drought + (1 | site)
##    Data: drought.prob2
## 
## REML criterion at convergence: 73.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3109 -0.6362 -0.3160  0.5727  1.9577 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.171    0.4135  
##  Residual             0.782    0.8843  
## Number of obs: 27, groups:  site, 11
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)   0.2778     0.3920 17.3623   0.709    0.488
## drought       0.4815     0.4519 19.1957   1.066    0.300
## 
## Correlation of Fixed Effects:
##         (Intr)
## drought -0.840
    legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Pop. size (after drought)", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ drought, drought.prob2, log = "y", main = "", las = 1)
    res <- lmer(log(larvae +1) ~ drought + (1|site), drought.prob2); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(larvae + 1) ~ drought + (1 | site)
##    Data: drought.prob2
## 
## REML criterion at convergence: 115.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.08430 -0.44864 -0.00434  0.34852  1.95495 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 2.770    1.664   
##  Residual             3.275    1.810   
## Number of obs: 27, groups:  site, 11
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)  
## (Intercept)    2.185      1.024 17.991   2.135   0.0468 *
## drought        1.785      1.137 22.835   1.570   0.1303  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## drought -0.798
    legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(juveniles+1) ~ drought, drought.prob2, log = "y", main = "", las = 1)
    res <- lmer(log(juveniles +1) ~ drought + (1|site), drought.prob2); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(juveniles + 1) ~ drought + (1 | site)
##    Data: drought.prob2
## 
## REML criterion at convergence: 99.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0094 -0.4241 -0.0583  0.4846  1.8962 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 2.006    1.416   
##  Residual             1.536    1.239   
## Number of obs: 27, groups:  site, 11
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)   0.5573     0.7783 19.2347   0.716    0.483  
## drought       1.9496     0.8404 24.4645   2.320    0.029 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## drought -0.773
    legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(frogs+1) ~ drought, drought.prob2, log = "y", main = "", las = 1)
    res <- lmer(log(frogs +1) ~ drought + (1|site), drought.prob2); summary(res)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(frogs + 1) ~ drought + (1 | site)
##    Data: drought.prob2
## 
## REML criterion at convergence: 75.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.97522 -0.41394 -0.00719  0.24147  2.13657 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.5121   0.7156  
##  Residual             0.6661   0.8161  
## Number of obs: 27, groups:  site, 11
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)   1.3769     0.4518 18.5096   3.048  0.00677 **
## drought       1.7862     0.5043 22.7391   3.542  0.00176 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## drought -0.802
    legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Dry versus not dry lake", side = 1, line = 2, outer = TRUE)